#include "opsimplex.h"

/*
 * simplex: resolve um pl pelo método simplex,
 * dada uma base inicial e os dados do problema. retorna se
 * encontrou o ótimo ou se o problema é ilimitado
 */
int simplex(PMatriz A, PMatriz c, PMatriz b, PMatriz x, PMatriz B, PMatriz invB, PMatriz cb, int* Ib, int* In, int m, int n, int* nit){
	if (A==NULL||c==NULL||x==NULL||B==NULL||invB==NULL||cb==NULL||Ib==NULL||In==NULL) return SIMPLEX_RET_ERR;
	if (m<=0||n<=0) return SIMPLEX_RET_ERR;
	if (A->nLin!=m||A->nCol!=n) return SIMPLEX_RET_ERR;
	if (c->nLin!=1||c->nCol!=n) return SIMPLEX_RET_ERR;
	if (b->nLin!=m||b->nCol!=1) return SIMPLEX_RET_ERR;
	if (x->nLin!=n||x->nCol!=1) return SIMPLEX_RET_ERR;
	if (B->nLin!=m||B->nCol!=m) return SIMPLEX_RET_ERR;
	if (invB->nLin!=m||invB->nCol!=m) return SIMPLEX_RET_ERR;
	if (cb->nLin!=1||cb->nCol!=m) return SIMPLEX_RET_ERR;


	if (verbose) fprintf(stderr,"starting simplex\n");

	PMatriz u, novaB, novaInvB, xBarraB;
	u = inicia(1,m);
	novaB = inicia(m,m);
	novaInvB = inicia(m,m);
	xBarraB = inicia(m,1);




	int cont = 0, i,k,ik,s,rc;
	double alphaK,cr;

	while (1){

		k = n;
		mult(cb,invB,u);
		mult(invB,b,xBarraB);

		if (verbose) fprintf(stderr,"updating variable values\n");
		for (i=0;i<m;i++)
			x->elems[Ib[i]][0]=xBarraB->elems[i][0]; //atualiza os valores das variáveis
		for (i=0;i<n-m;i++)
			x->elems[In[i]][0]=0.0;

		for (i=0; i<n-m; i++){
			cr = calculaCustoReduzido(In[i],c,u,A);
			if (!equals(cr,0.0) && cr<0 && In[i]<k){ //regra de bland
				k = In[i];							 //pegue o de menor indice
				ik = i;
			}
		}
		if (verbose) fprintf(stderr,"selected k=%d\n",k);
		if (k == n) {rc = SIMPLEX_RET_OT; break;} //não há variável em In com custo negativo
		cont++;
		if (verbose) fprintf(stderr,"iteration %d\n",cont);

		s = calculaAlphaK(k,invB,A,xBarraB,Ib,&alphaK);
		if (verbose) fprintf(stderr,"selected s=%d: alphak=%lf\n",s,alphaK);
		if (s==-1) {rc = SIMPLEX_RET_IL; break;} //não há limite superior para k






		trocaCol(B,s,A,k,novaB);
		inv(novaB,B,invB,s,novaInvB); //atualiza as matrizes B
		freeMatriz(B);                //e invB
		B = inicia(m,m);
		clone(novaB,B);
		freeMatriz(invB);
		invB = inicia(m,m);
		clone (novaInvB,invB);
		novaB = NULL;
		novaInvB = NULL;
		novaB = inicia(m,m);
		novaInvB = inicia(m,m);

		if (verbose) fprintf(stderr,"base updated\n");


		cb->elems[0][s] = c->elems[0][k]; //atualiza cb
		int aux = Ib[s]; //atualiza Ib e In
		Ib[s] = k;
		In[ik] = aux;



	}

	freeMatriz(u);
	freeMatriz(novaB);
	freeMatriz(novaInvB);
	freeMatriz(xBarraB);

	*nit = cont;
	return rc;
}

/*
 * simplexDuasFazes: resolve um pl pelo método simplex de duas
 * fases, a partir dos dados do problema. retorna se encontrou o
 * ótimo ou se o problema é ilimitado ou inviável
 */
int simplexDuasFazes(PMatriz A, PMatriz c, PMatriz b, PMatriz x){
	if (A==NULL||c==NULL||b==NULL) return SIMPLEX_RET_ERR;
	int m=A->nLin, n=A->nCol;
	if (c->nLin!=1||c->nCol!=n) return SIMPLEX_RET_ERR;
	if (b->nLin!=m||b->nCol!=1) return SIMPLEX_RET_ERR;


	if (verbose) fprintf(stderr,"starting auxiliary structures for first phase\n");

	it1 = 0;
	it2 = 0;

	time_t tini, tfim;
	int i,rc,rci,j,aux, nit;
	PMatriz B, invB, artA, artC, artX, cb, temp1, temp2, temp3, sol, retInvB;
	B = inicia(m,m);
	invB = inicia(m,m);
	artA = inicia(m,m+n);
	artC = inicia(1,m+n);
	artX = inicia(m+n,1);
	freeMatriz(x);
	x = inicia(n,1);
	cb = inicia(1,m);
	temp1 = inicia (1,m);
	temp2 = inicia(m,m);
	temp3 = inicia(m,m);
	sol = inicia(1,1);

	if (verbose) fprintf(stderr,"checking for negative resources\n");
	for (i=0;i<m;i++)
		if (b->elems[i][0]<0.0 && !equals(b->elems[i][0],0.0)){
			b->elems[i][0]*=-1;
			for (j=0;j<n;j++)
				A->elems[i][j]*=-1;
		}



	newId(B,m);
	newId(invB,m);
	append(A,B,artA);

	for(i=0;i<n;i++)
		artC->elems[0][i] = 0;
	for(i=0;i<m;i++){
		artC->elems[0][i+n] = -1;
		cb->elems[0][i] = -1;
	}

	for(i=0;i<n;i++)
		artX->elems[i][0] = 0;
	for(i=0;i<m;i++)
		artX->elems[i+n][0] = b->elems[i][0];

	int *Ib, *artIn, *In;
	(Ib = (int*) malloc(m*sizeof(int))) || die(MSG_ERR_MALLOC);
	(artIn = (int*) malloc(n*sizeof(int))) || die(MSG_ERR_MALLOC);
	(In = (int*) malloc((n-m)*sizeof(int))) || die(MSG_ERR_MALLOC);

	for (i=0;i<m;i++)
		Ib[i]=n+i;
	for (i=0;i<n;i++)
		artIn[i]=i;

	if (verbose) fprintf(stderr,"starting first phase\n");
	tini = time(NULL);

	rci = simplex(artA,artC,b,artX,B,invB,cb,Ib,artIn,m,m+n,&nit); //primeira fase
	it1 = nit;

	tfim = time(NULL);
	t1 = difftime(tini,tfim);
	if (verbose) fprintf(stderr,"first phase done: rc=%d\n",rci);


	if (rci!=SIMPLEX_RET_OT) rc = rci;
	else {
		mult(cb,invB,temp1);
		mult(temp1,b,sol);
		if(!equals(getEscalar(sol),0.0)) rc = SIMPLEX_RET_INV;


		else{

			if (verbose) fprintf(stderr,"removing non basic variables\n");

			j=-1;
			for(i=0;i<m;i++)
				if(Ib[i]>=n){                            //troca possiveis variaveis
					while (j<n-1 && artIn[++j]>=n);      //não básicas na base
					if (artIn[j]>=n) die(MSG_ERR_BASE);  //devido a soluções degeneradas
					aux = artIn[j];
					artIn[j] = Ib[i];
					Ib[i] = aux;
					trocaCol(B,i,A,artIn[j],temp2);
					inv(temp2,B,invB,i,temp3);
					freeMatriz(B);
					freeMatriz(invB);
					B=temp2;
					invB=temp3;
				}

			for (i=0;i<n-m;i++){
				while (j<n-1 && artIn[++j]>=n);
				if (artIn[j]>=n) die(MSG_ERR_BASE);
				In[i]=artIn[j];
			}


			if (verbose) fprintf(stderr,"updating costs and variable values\n");

			for (i=0;i<n;i++)
				x->elems[i][0] = artX->elems[i][0];

			for (i=0;i<m;i++)
				cb->elems[0][i] = c->elems[0][Ib[i]];


			if (verbose) fprintf(stderr,"starting second phase\n");
			tini = time(NULL);

			rc = simplex(A,c,b,x,B,invB,cb,Ib,In,m,n,&nit);
			it2 = nit;

			tfim = time(NULL);
			t2 = difftime(tini,tfim);
			if (verbose) fprintf(stderr,"second phase done: rc=%d\n",rc);
		}
	}

	freeMatriz(B);
	freeMatriz(invB);
	freeMatriz(artA);
	freeMatriz(artC);
	freeMatriz(artX);
	freeMatriz(cb);
	freeMatriz(temp1);
	freeMatriz(temp2);
	freeMatriz(temp3);
	freeMatriz(sol);

	if (Ib) free(Ib);
	if (artIn) free(artIn);
	if (In) free(In);

	return rc;


}

/*
 * calculaAlphaK: calcula o alfaK. retorna s, a variável que sai
 * da base, e o valor de alfaK vai por referência
 */
int calculaAlphaK(int pk, PMatriz invB, PMatriz A, PMatriz xBarraB, int*Ib, double* alphaK){
	if(invB==NULL||A==NULL||xBarraB==NULL||Ib==NULL||pk>=A->nCol) return -2;
	int m = invB->nCol;
	if(!isQuadrada(invB)||(xBarraB->nCol!=1)||xBarraB->nLin!=invB->nLin||A->nLin!=invB->nLin) return -2;

	int s,i;
	double cand;
	PMatriz yk, aj;
	yk = inicia(m,1);
	aj = inicia(m,1);

	getColuna(pk,A,aj);
	mult(invB,aj,yk);

	s=-1;
	*alphaK = HUGE_VAL;
	for (i=0;i<invB->nCol;i++){
		if(yk->elems[i][0]>0.0 && !equals(yk->elems[i][0],0.0)){
			cand = xBarraB->elems[i][0]/yk->elems[i][0];
			if(cand<*alphaK && !equals(cand,*alphaK)){ //estritamente menor
				*alphaK = cand;                       //regra de bland
				s = i;                               //pega o menor i
			}
		}
	}

	freeMatriz(yk);
	freeMatriz(aj);
	return s;
}

/*
 * calculaCustoReduzido: calcula o custo reduzido associado a
 * variável j
 */
double calculaCustoReduzido(int pj, PMatriz c, PMatriz u, PMatriz A){
	if (c==NULL||u==NULL||A==NULL||pj>=A->nCol) return (-HUGE_VAL);
	PMatriz aj, uaj;
	int m = A->nLin;
	aj = inicia(m,1);
	uaj = inicia(1,1);

	if  (
		(!getColuna(pj,A,aj))||
		(!mult(u,aj,uaj))||
		(!isEscalar(uaj))
		) return (-HUGE_VAL);

	double result = (getEscalar(uaj))-(c->elems[0][pj]);

	freeMatriz(aj);
	freeMatriz(uaj);
	return result;
}


